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1 Abstract 

A glimpse is provided of the research program in 
stability, transition and turbulence based on numer- 
ical simulations. This program includes both the so- 
called abrupt and the restrained transition processes. 
Attention is confined to the prototype problems of 
channel flow and the parallel boundary layer in the 
former category and the Taylor-Couette flow in the 
latter category. It covers both incompressible flows 
and supersonic flows. Some representative results are 
presented. 


3 Introduction 

The phenomena of transition to turbulence are so 
complex as to defy a unified theory at the present 
time. As such, direct numerical simulations of these 
phenomena within the framework of Navier-Stokes 
equations have assumed a dimension equal in im- 
portance to experimental and even theoretical stud- 
ies. An effective simulation not only mimics a phys- 
ical experiment but has the added advantages of of- 
fering readily retrievable clean information wherever 
and whenever it is needed, and also a precise control 
of parameters unachievable in the physical experi- 
ment. Theoretical and mathematical approaches pro- 
vide progress in understanding through the processes 
of abstraction and idealization; the results of anal- 
ysis then furnish, apart from specific predictions, a 
deeper comprehension of the underlying general prin- 
ciples. These features are not intrinsic to a numerical 
simulation. However, they can be built into the simu- 
lations to some extent. Careful considerations of such 
factors are the bases of the research program in this 
area. Here at LaRC/ICASE our research program in- 
cludes both the so-called abrupt and restrained tran- 
sition processes. Attention is confined to the proto- 
type problems of the channel flow and the parallel 
boundary layer in the former category, and the wake 
and the Taylor-Couette flows in the latter category. 
It covers both incompressible flows and supersonic 
flows. 
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2 Nomenclature 

C p : specific heat at constant pressure 

C v : specific heat at constant volume 

D : average of diagonal stress components 

F : forcing function 

M : Mach number 

OD : average of off-diagonal stress components 
p : pressure 

Pr : Prandtl number 

Re : Reynolds number 

t : time 

T : temperature 

u : fluid velocity 

7 : ratio of specific heats (C p /C v ) 

6 : Kronecker tensor 

V : gradient operator 

k : thermal conductivity 

fx : molecular viscosity 

$ : viscous dissipation function 

p : density 

S : average of stress components 

on the scalar level 
<r : viscous stress tensor 

9 : (T-T^HT-T^) 

V : average of stress components 

on the vector level 


Superscripts 

: Favre average 

Subscripts 

k , l : Cartesian indices 

w : wall value 

oo : freestream value 



4 Poiseuille and Blasius Flows 

4.1 Background 

Although our knowledge of laminar-turbulent tran- 
sition is by no means complete, the broad features 
are now clear at least in low subsonic shear flows. 
This is in a large measure due to the the clas- 
sical experiments of Klebanoff, Tidstrom and Sar- 
gent [l], Kovasznay, Komoda and Vasudeva [2], Hama 
and Nutant [3], These experiments were conducted 
in a controlled, identifiable disturbance environment 
where the nondeterministic disturbances were im- 
measurably small. Briefly stated, the instability and 
transition process involves the following stages: 1) 
a primary instability evolving in accordance with 
the linear theory, 2) a secondary instability lead- 
ing to the emergence of a flow of a pronounced 
three-dimensional nature, and the appearance of the 
streamwise vortex system, 3) the development of de- 
tached high-shear layers, 4) at least a tertiary in- 
stability (if not an extremely rapid succession of a 
sequence of instabilities) resulting in turbulent spot 
development and 5) convection and coalescence of 
turbulent spots to form a fully developed turbulent 
flow. Similar patterns have been observed in plane 
Poiseuille flow by Nishioka, Iida and Ichikawa [4]. 
More recently, the experiments of Kachanov and 
Levchenko [5], and Saric and Thomas [6] have un- 
covered some new details in the secondary instability 
stage of the transition process in boundary layers. 
The flow pattern outlined above is relevant to natu- 
ral transition under some limited, but not all, condi- 
tions. The reviews of Stuart [7] and Tani [8] provide 
a detailed discussion of the classical experiments and 
the then available theoretical explanations of the iso- 
lated stages or events of the transition process. The 
more recent review of Herbert and Morkovin [9] pro- 
vides an overview of applied problems in stability and 
transition, and a discussion of the various theoreti- 
cal approaches to the secondary instability problem. 
Craik’s monograph [10] deals with the general prob- 
lem of wave interactions in various fields, and includes 
a discussion of the secondary instability pertinent to 
wall-bounded shear flows. 

The primary instability has as its basis linear the- 
ory and has long been well established. The the- 
ory for the secondary instability is however more re- 
cent, but is also fairly well established (Herbert [11], 
Nayfeh [12]). There have been no improvements in 
the theoretical models for the other stages since the 
survey article of Stuart [7] Full scale numerical sim- 
ulations provide a unified basis for investigating the 


relevance and relative importance of various events 
right up to the formation of turbulent spots. In the 
present paper we shall dwell on numerical experi- 
ments which have uncovered the secondary instability 
associated with the center modes in plane Poiseuille 
flow. We shall also discuss simulations which have 
brought out the relative effectiveness of the various 
laminar flow control (LFC) techniques such as heat- 
ing, favorable pressure gradient and suction. We will 
present some representative results in each category. 

The incompressible variable property Navier- 
Stokes equations, in usual notation, are (Zang and 
Hussaini [13]) 

V * u = 0 (1) 

+ u Vu = -Vp+ i_V • ( M Vu) + Fu (2) 

§r + " < 3 > 

where the lengths are scaled by the half channel- 
width (or for the boundary layer by the displacement 
thickness at the streamwise station of interest), ve- 
locities by the corresponding center velocity for the 
channel flow (or the corresponding mainstream ve- 
locity for the boundary layer), and pressure by the 
dynamic head; the forcing functions F \ i and F$ are 
designed to ensure parallel flow. The viscosity and 
conductivity are given functions of temperature for 
the heated wall case, and constants in other cases. 

The applied boundary conditions assume periodic- 
ity in the horizontal directions, and the no-slip ve- 
locity on solid walls. For the heated water bound- 
ary layer case, the no-slip condition for the velocity 
and a uniform temperature are enforced at the wall; 
perturbations are assumed to be zero in the main- 
stream. This permits a Fourier discretization in the 
streamwise and spanwise directions, and a Cheby- 
shev discretization in the vertical direction. The ini- 
tial condition consists of a triad of waves - one two- 
dimensional wave and two skewed waves. The general 
fractional step algorithm for efficiently solving the rel- 
evant equations is given by Zang and Hussaini [14]. 
A survey of spectral algorithms for fluid dynamic 
calculations is given in Hussaini and Zang [15], and 
Canuto, Hussaini, Quarteroni and Zang [16]. The lat- 
ter is a monograph which deals with both theory and 
applications of spectral methods. 

4.2 Instability Due To A Triad Of 
Center Modes 

The so-called center modes of plane channel flows are 
in fact temporal eigenfunctions of the Orr-Sommer- 
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feld equation which decay in time. However, a com- 
bination of a two-dimensional center mode with two 
skewed modes is susceptible to an instability simi- 
lar to the secondary instability encountered in a sim- 
ilar non-resonant triad of wall modes. The essen- 
tial qualitative difference between the center mode 
and the wall mode (as is evident from Fig. la and 
lb) is that the maximum amplitude for the former 
occurs away from wall towards the channel center, 
while for the latter it occurs near the wall, and hence 
the terminology. There are of course other quanti- 
tative differences with regard to wavelengths, phase 
velocities and decay rates. The Reynolds number 
for the simulation was 8000 based on the channel 
half width, and both the streamwise and spanwise 
wavenumbers were unity; the initial two-dimensional 
and three-dimensional amplitudes were 10% and 3% 
respectively of the channel center velocity. The reso- 
lution requirements were monitored and the grid was 
refined to resolve the fine structures as they emerged. 
The finest grid was 96x128x192. Plotted in Fig. 2 is 
the time history of the harmonic contents of the so- 
lution. Note that the two-dimensional mode, labeled 
(1,0), decays almost exponentially, but the three- 
dimensional mode, labeled (1,1), grows after an initial 
decay for about three periods 2 it being the charac- 
teristic period. The results are presented at t = 87 
(a little less than 14 periods) at which time the finest 
grid was just about sufficient to resolve the sharp gra- 
dients of the flow field. The vortex lines are displayed 
in Fig. 3a, and for comparison purposes vortex lines 
for an analogous case of wall modes are presented 
in Fig. 3b. The streamwise and spanwise vorticity 
contours are displayed in Fig. 4 at four equidistant 
streamwise planes over a wavelength. Fig. 5 dis- 
plays similar contours for the case of the wall modes. 
The qualitative similarities and differences are obvi- 
ous. The center-mode vortex loop is comparatively 
away from the wall, and the pinching at the neck is 
less intense than for the wall-mode vortex loop. Also, 
the high-shear layer associated with the center modes 
is less intense than the one associated with the wall 
modes. The physical relevance of the center mode 
instability to the high-intensity bypass to transition 
is as yet undetermined. 

4.3 Effect of Heating on the Second- 
ary Instability of Blasius Flow 

Laminar flow control (LFC) vehicles might operate 
in a finite-amplitude disturbance environment, so it 
is of interest to examine the relative performance of 
the various LFC techniques in the nonlinear regime. 


While some experimental studies have investigated 
the influence of nonlinear disturbance sources such as 
roughness, vibrations, free-stream turbulence, etc. on 
the basic transition process, there have been no com- 
parative studies on the sensitivity of the LFC tech- 
niques such as heating, pressure gradient and suction 
in the nonlinear regime. In the absence of any rele- 
vant nonlinear stability theories due to the complex- 
ity of the problem, numerical simulations can play 
a crucial role in this field. As an example of such a 
use of simulation, we have performed a Navier-Stokes 
calculation involving conditions similar to those used 
in the experiment of Kovasznay, Komoda and Va- 
sudeva [2] except that in our simulation the wall was 
heated. The idealization of a parallel boundary layer 
was used mainly in order to meet the resolution re- 
quirements while keeping the computing time within 
reasonable limits. The Reynolds number of the simu- 
lation was 1100 based on the displacement thickness, 
and the initial amplitudes of the two-dimensional and 
three-dimensional Tollmien-Schlichting waves (wall 
modes) were respectively 2.7% and 0.8% of the free 
stream velocity. The finest grid was 72x162x192. 

Three different situations were studied: 1) the un- 
controlled case, 2) the heated fixed temperature case, 
and 3) the heated active temperature case. In the 
heated fixed temperature case, the temperature was 
kept fixed at the initial value pertinent to the mean 
flow conditions, and the temperature evolution was 
totally neglected. In the heated active temperature 
case, the temperature evolution was taken into ac- 
count by solving the temperature equation along with 
the momentum equations. In both the cases the wall 
temperature was 2.75% above the free stream tem- 
perature. Fig. 6 shows the harmonic history of 
the perturbation energy. Note that the (1,0) mode 
which grows in the uncontrolled case decays in the 
heated almost up to 4 periods. Heating damps the 
(1,1) mode initially, but it starts growing within the 
first period. It appears that when the energy in the 
three-dimensional wave overtakes that in the two- 
dimensional wave, it feeds energy partly into the two- 
dimensional wave. Fig. 7 shows the spanwise vortic- 
ity contours on the peak plane. In the uncontrolled 
case (Fig. 7 top left) a kink develops in the high-shear 
layer at time t equal to three Tollmien-Schlichting 
periods. It is generally accepted that a irrevocably 
quick succession of events follows thereafter leading 
to a turbulent spot formation. Heating the wall to 
2.75% above the free stream temperature diffuses the 
high-shear layer as is obvious from Fig. 7 (bottom 
left). However, within the subsequent one and one 
fourth period, turbulent spot formation appears to 
become imminent (Fig. 7 top right). In the fixed 
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temperature case, it is clear from Fig. 7 (bottom 
right) that the high-shear layer formation is mellowed 
down even up to four and one fourth periods. The 
fixed temperature case overpredicts the weakening ef- 
fect of heating on the secondary instability. In other 
words, the effect of temperature evolution is signifi- 
cant and deleterious in the nonlinear regime whereas 
it is quite negligible in the linear regime. 


5 Taylor-Couette Flow 

The instability and abrupt transition process are typ- 
ical of aerodynamic flows. On the other hand, the 
transition to turbulence by spectral evolution are typ- 
ical of geophysical fluid dynamics where there is a 
destabilizing force field. The Taylor-Couette flow 
in the annulus of concentric rotating circular cylin- 
ders typifies a class of geophysical problems involv- 
ing instabilities and turbulence due to the presence 
of a destabilizing centrifugal force field. DiPrima and 
Swinney [17] provide an excellent review of instabil- 
ities and transitions in Taylor-Couette flows. Our 
interest lies in finite-length Taylor-Couette flow as 
it offers the situation of weak turbulence which is 
amenable to direct simulation. The prime objective 
of our research effort to examine the various quali- 
tative mathematical theories (Benjamin [18]). Some 
of the routes to mathematical chaos appears to have 
been observed in these flows. The chaos theories are 
based on model equations. Our investigations will 
establish the relevance of these theories to Navier- 
Stokes equations. 

This problem, unlike those considered in the pre- 
ceding sections, is truly inhomogeneous in two di- 
rections. A spectral algorithm necessarily involves 
Chebyshev polynomial expansions in radial and axial 
directions, and Fourier expansion in the azimuthal di- 
rection. The fractional step algorithm used to study 
the channel flows and the parallel boundary layers 
is easily extended to treat this problem (Streett and 
Hussaini [19]). 

The initial phase of study focussed on the steady- 
state bifurcations in axisymmetric Taylor-Couette 
flow (with the inner cylinder rotating, and the outer 
cylinder and end walls remaining stationary) for as- 
pect ratios of order unity. These problems provide 
a stringent test for a time-accurate method designed 
to simulate the delicate unsteady processes leading 
to transition. Among the numerous simulations car- 
ried out to compare with the available experiments, 
we present some typical results for the geometry of 
Benjamin and Mullin [20] which had an aspect ratio 


of 1.05 and radius ratio of 0.615. The symmetric two- 
cell mode is established at a relatively low Reynolds 
number of 62, and then the Reynolds number is im- 
pulsively raised to 150. The development of the asym- 
metric single-cell flow represented by the streamlines 
in a cross-sectional plane is shown in Fig. 8, and Fig. 
9 displays the order parameter (which is a measure 
of the asymmetry of the flow) and energy as func- 
tions of time in units of characteristic diffusion time 
scale. The continuing simulations are on the verge 
of capturing bifurcations into time-dependent states. 
The final objective is to calculate the dimension of 
the strange attractor which is supposed to represent 
the weak turbulence. 


6 Compressible Transition and 
Turbulence 


Transition to turbulence in supersonic and hyper- 
sonic flows is a gray area. There are differing experi- 
mental results, and unexplained visual observations. 
The picture is quite piecemeal compared to low Mach 
number flows. The purpose of the research program 
is to answer questions such as: a) do secondary in- 
stability mechanisms observed in low subsonic flows 
persist in supersonic flows, b) is the turbulent spot 
formation the usual way transition occurs, and c) is 
the transition process abrupt or restrained. The ini- 
tial focus of our research program is on question (a), 
i.e., whether the known prototypes of instabilities in 
incompressible flows persist at high Mach numbers. 
A spectral simulation of a three-dimensional insta- 
bility process in a parallel boundary layer at Mach 
number 4.5 was carried out. 

The full compressible Navier-Stokes equations are, 
in dimensionless form, (Erlebacher and Hussaini [21]) 
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and the equation of state is 

7 M£,p = pT (9) 

with 7 = C p /C v , the ratio of specific heats. At 
high Mach number, temperature variations across 
the boundary layer become important and the tem- 
perature dependence of viscosity must be taken into 
account. The Prandtl number was assumed equal 
to 0.7, and the Sutherland’s law was prescribed for 
viscosity variation with the temperature. The free 
stream Mach number was 4.5, and the Reynolds num- 
ber was 10000 based on the displacement thickness. 
The boundary conditions were identical to those used 
in the incompressible boundary layer simulations dis- 
cussed in the earlier section. The initial condition 
consisted of an unstable two-dimensional wave and a 
pair of three-dimensional oblique waves superimposed 
on a prescribed parallel mean flow. These waves were 
the eigenfunctions of an eighth order eigenvalue prob- 
lem, the compressible analogue of the incompressible 
Orr-Sommerfeld equation. The initial amplitudes of 
the two-dimensional and three-dimensional perturba- 
tions were taken to be respectively 5.4% and 1.2% 
of the free stream velocity. The temporal evolu- 
tion of this triad of waves was followed for about 
eight periods of the initial two-dimensional Tollmien- 
Schlichting wave (TS wave). The time history of sev- 
eral Fourier harmonics of the perturbation kinetic en- 
ergy (integrated in the vertical direction) is shown 
in Fig. 10. After a period of slow growth of the 
three-dimensional wave, labeled (1,1) mode, the non- 
linear interactions between the (0,1) spanwise mode 
and the primary (1,0) mode trigger the secondary in- 
stability after approximately 5 TS periods. Contour 
plots of the streamwise velocity (Fig. 11) at the time 
of 5 TS periods illustrate the incipient stages of a K- 
type breakdown (Erlebacher and Hussaini [21]). In 
the supersonic simulations, the critical layer is an or- 
der of magnitude further from the wall than in the 
incompressible cases. One consequence is that the 
structures are farther removed from the walls. Ini- 
tial results indicate that compressibility retards the 
onset of the secondary instability. However, a more 
extensive parameter study is required before any firm 
conclusion can be drawn. 

Our research program has recently focused its at- 
tention on the development of subgrid-scale mod- 
els for compressible turbulence. A compressible ex- 
tension of Bardina’s linear combination model [22] 
has been developed. Favre- averaging was employed 
to produce spatially averaged compressible Navier- 
Stokes equations which closely resemble the incom- 
pressible ones. For compressible as opposed to in- 
compressible flow, the trace of the subgrid stresses 


cannot be incorporated with the pressure since the 
pressure is a true thermodynamic variable. These 
isotropic stresses must now be modeled, and a new 
model constant, termed the isotropic constant, must 
be included. 

This new model has been tested against direct sim- 
ulations at Mach numbers ranging from 0.0 to 0.6 and 
on grids of 64 3 to 128 3 . The turbulence model was 
compared against this direct simulation (DS) data 
suitably injected onto a coarser grid. The coarse 
grid would normally be used in a large eddy simula- 
tion. A least square fit between the total stress calcu- 
lated from the DS and the modeled stress computed 
from spatially filtered velocities allows the constants 
to be determined. The constant determination was 
performed on the vector level (Erlebacher, Hussaini, 
Speziale, Zang [23]). Table 1 summarizes the correla- 
tion coefficients between the modeled and the exact 
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Table 1: 


total stresses at Mach 0 and Mach 0.6. Results are 
presented on the tensor, vector and scalar levels. The 
vector level is the divergence of the tensor, while the 
scalar level is defined by the product of V * u and u. 
The Leonard stress is omitted from the total subgrid 
stress because it can be calculated exactly. Corre- 
lation coefficients of over 80% are obtained on the 
tensor level, and above 70% on both the vector and 
scalar levels. A more thorough study can be found in 
Erlebacher et al. [23]. 

7 Concluding Remarks 

Numerical algorithms for studying the physics of 
transition and turbulence in simple geometries have 
been developed. Detailed studies of the highly nonlin- 
ear stages of transition prior to turbulent spot forma- 
tion have been made. The use of simulations to study 
the sensitivity of laminar flow control techniques in 
the nonlinear regime has been demonstrated. A 
new instability mechanism associated with the center 
modes in channel flows has been uncovered, although 
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its physical significance is as yet undetermined. The 
first numerical simulations of their kind are being car- 
ried out for supersonic boundary layer flows. The pre- 
liminary results revealed the existence of a secondary 
instability similar to the one in incompressible flows. 
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Figure 4. Streamwise (left) and spanwise (right) 
vorticity contours at t = 87 for a channel center 
mode at a Reynolds number of 8000. The contour 
interval is 0.20. Dashed lines indicate negative 
contours. 



Figure 5. Streamwise (left) and spanwise (right) 
vorticity contours at t = 182 for a channel wall 
mode at a Reynolds number of 8000. The contour 
interval is 0.40. Dashed lines indicate negative 
contours. 



Figure 6. Harmonic energy history of an 
uncontrolled (solid lines) and a heated (dashed 
lines) boundary layer undergoing a K-type 
transition at a Reynolds number of 1100. 
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Figure 8. Cross-sectional streamline at 5 stages Figure 10. Harmonic energy history for a Mach 4.5 
of two-cell/one-cell exchange process. boundary layer undergoing transition at a Reynolds 

number of 10,000. 
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Figure 11. Streamwise vorticity contours for the 
compressible boundary layer transition at Mach 
4.5. The contour interval is 0.20. Dashed lines 
indicate negative contours. 
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Figure 9. Time history of order parameter and 
perturbation energy. 
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